************************************************** 
*
*  Post-estimation analysis
*  
*  This file
*
*  1. conducts specification tests (LR-tests)
*     as stated in the paper and SI 
*  2. produces fit statistics (information 
*     criteria) for models 1-3 as stated in the
*     paper and SI
*  3. produces Table SI-2
*  4. produces Figure 2
*
**************************************************



* Resurrect results matrices

cd "..."

#delimit ;

local varsM1 "dist_med_part dist_med_kpart d_alter ant_eigReg ant_eigKonf 
	ant_eigStand ant_eigBeruf ant_eigVorp ant_eigParlAmt ant_eigPolAmt
	ant_eigVerein ant_eigStrafe" ;
local varsM2 "dist_med_part dist_med_kpart" ;
local varsM3 "d_alter ant_eigReg ant_eigKonf ant_eigStand ant_eigBeruf 
	ant_eigVorp ant_eigParlAmt ant_eigPolAmt ant_eigVerein ant_eigStrafe" ;

#delimit cr

local clubs "Do De Wh Wü Au La Ca Mi" 

forvalues i = 1/3 {
	use estimates_M`i'.dta, clear
	mkmat _all, matrix(M`i')
	matrix rownames M`i' = `varsM`i'' `clubs'
}

* Load data and results
use data_and_results.dta, clear

* Load model 3 results
est use model3




* Specification tests reported in paper
*****************************************

* LR-test: b_part = b_kpart (for model 1 and model 2)

forvalues i = 1/2 {

	quietly sum logL_M`i'
	scalar ll_M`i' = r(mean)

	quietly sum logL0_M`i'
	scalar ll0_M`i' = r(mean)

	display as text "LR test of beta:dist_med_part " /// 
		"= beta:dist_med_kpart in Model " `i' as text ":" _newline ///
		"p > chi2(1, 2*(" ll_M`i' as text "-" ll0_M`i' ///
		as text ")) = " as result 1 - chi2(1, 2 * (ll_M`i' - ll0_M`i'))
}

* LR-test: model 1 = model 2

display as text "LR test of M1 = M2:" _newline ///
	"p > chi2(10, 2*(" ll_M1 as text "-" ll_M2 ///
	as text")) = " as result 1 - chi2(10, 2 * (ll_M1 - ll_M2))

* LR-test: model 1 = model 3

est restore mod3
scalar ll_M3 = e(ll)

display as text "LR test of M1 = M3:" _newline ///
	"p > chi2(2, 2*(" ll_M1 as text "-" ll_M3 ///
	as text")) = " as result 1 - chi2(2, 2 * (ll_M1 - ll_M3))



	
* Measures of model fit
*************************

* Formulas for AIC and BIC were taken from:

help bic_note

* To calculate BIC, we use the number of MPs (N = 506). 
* Stata's default is to use the number of observations (N = 4554).

est restore mod3
scalar ll_M3 = e(ll)

* AIC and BIC

matrix k = (20\ 10\ 18)			// number of model parameters 

forvalues i = 1/3 {
	scalar aic_M`i' = -2 * ll_M`i' + 2 * k[`i', 1]
	scalar bic_M`i' = -2 * ll_M`i' + log(506) * k[`i', 1]	
}




* Table SI-2: Regression results  
**********************************

* Matrix to fill with regression table
matrix etable = J(43, 3, .)

* Model 1

* Coefs 
local m = 0
forvalues n = 1(2)39 {
	local ++m
	matrix etable[`n', 1] = M1[`m', 1]
}

* SEs
local m = 0
forvalues n = 2(2)40 {
	local ++m
	matrix etable[`n', 1] = M1[`m', 2]
}

* Model 2

* Coefs and SEs for ideology 
matrix etable[1, 2] = M2[1, 1]
matrix etable[2, 2] = M2[1, 2]
matrix etable[3, 2] = M2[2, 1]
matrix etable[4, 2] = M2[2, 2]

* Constants incl SEs 
local m = 2
forvalues n = 25(2)39 {
	local ++m
	matrix etable[`n', 2] = M2[`m', 1]
}
local m = 2
forvalues n = 26(2)40 {
	local ++m
	matrix etable[`n', 2] = M2[`m', 2]
}

* Model 3

* Coefs  
local m = 0
forvalues n = 5(2)39 {
	local ++m
	matrix etable[`n', 3] = M3[`m', 1]
}

* SEs
local m = 0
forvalues n = 6(2)40 {
	local ++m
	matrix etable[`n', 3] = M3[`m', 2]
}

* N, AIC, BIC
matrix etable[41, 1] = 4554, 4554, 4554
matrix etable[42, 1] = aic_M1, aic_M2, aic_M3
matrix etable[43, 1] = bic_M1, bic_M2, bic_M3

* exporting results (numbers only) to be edited with word processor
esttab matrix(etable, fmt(3)) using regression_table.rtf, replace

 
 
	
* Figure 2: Coefficient plot 
******************************

* b/w
# delimit ;
coefplot 
	(matrix(M1[,1]), ci((M1[,3] M1[,4])) label("Full model")) 						
	(matrix(M2[,1]), ci((M2[,3] M2[,4])) label("Ideological similarity only"))		
	(matrix(M3[,1]), ci((M3[,3] M3[,4])) label("Without ideological similarity") 
		xline(0) legend(rows(3) span ring(0) pos(5)) 								
		note("Conditional logit coefficients with 95 percent confidence intervals." 
		"Coefficients are directly comparable because all predictors have a range of 1.", span)
		scheme(s1mono))	
	, ylabel(, noticks)	xlabel(-5(5)20)	xscale(range(23))
		coeflabel(dist_med_part ="Ideology (alternatives: parties)"
		dist_med_kpart ="Ideology (alternative: unaffiliated)" 	
		Do ="Donnersberg" De = "Deutscher Hof"						
		Wh = "Westendhall" Wü = "Württemberger Hof" 				
		Au = "Augsburger Hof" La = "Landsberg"
		Ca = "Casino" Mi = "Milani") 			
		order(dist_med_part dist_med_kpart d_alter ant_eigReg
		ant_eigKonf ant_eigStand ant_eigBeruf ant_eigVorp ant_eigParlAmt
		ant_eigPolAmt ant_eigVerein ant_eigStrafe) 						
		headings(dist_med_part = "{bf: Similar Ideology}" 				
		d_alter = "{bf: Similar Cultural Background}" 					
		ant_eigStand = "{bf: Similar Socio-Economic Status}" 			
		ant_eigVorp = "{bf: Similar Political Experience}" 				
		Do = "{bf: Alternative-specific constants}", gap(.5)) scale(*.7) ;
#delimit cr
graph export coefplot.emf, replace


* color
# delimit ;
coefplot 
	(matrix(M1[,1]), ci((M1[,3] M1[,4])) label("Full model")) 						
	(matrix(M2[,1]), ci((M2[,3] M2[,4])) label("Ideological similarity only"))		
	(matrix(M3[,1]), ci((M3[,3] M3[,4])) label("Without ideological similarity") 
		xline(0) legend(rows(3) span ring(0) pos(5)) 								
		note("Conditional logit coefficients with 95 percent confidence intervals." 
		"Coefficients are directly comparable because all predictors have a range of 1.", span)
		scheme(s1color))	
	, ylabel(, noticks)	xlabel(-5(5)20)	xscale(range(23))
		coeflabel(dist_med_part ="Ideology (alternatives: parties)"
		dist_med_kpart ="Ideology (alternative: unaffiliated)" 	
		Do ="Donnersberg" De = "Deutscher Hof"						
		Wh = "Westendhall" Wü = "Württemberger Hof" 				
		Au = "Augsburger Hof" La = "Landsberg"
		Ca = "Casino" Mi = "Milani") 			
		order(dist_med_part dist_med_kpart d_alter ant_eigReg
		ant_eigKonf ant_eigStand ant_eigBeruf ant_eigVorp ant_eigParlAmt
		ant_eigPolAmt ant_eigVerein ant_eigStrafe) 						
		headings(dist_med_part = "{bf: Similar Ideology}" 				
		d_alter = "{bf: Similar Cultural Background}" 					
		ant_eigStand = "{bf: Similar Socio-Economic Status}" 			
		ant_eigVorp = "{bf: Similar Political Experience}" 				
		Do = "{bf: Alternative-specific constants}", gap(.5)) scale(*.7) ;
#delimit cr
graph export coefplot_color.emf, replace

